library(plotly)
library(readr)
library(gganimate)
library(ggmap)
library(dplyr)
library(ggplot2)
library(highcharter)
library(ISOcodes)
library(countrycode)
library(stringr)
library(lubridate)


worldwide = read_csv("data_eq/worldwide.csv")
worldwide$year = format(worldwide$time, "%Y")
for (i in 1:nrow(worldwide)) {
  worldwide$place[i] = worldwide$place[i] %>% str_split(", ") %>% .[[1]] %>% tail(1)
}
historical = read_rds("data_eq/historical_web_data_26112015.rds")
disaster = read_delim("data_eq/disaster.txt", "\t", escape_double = FALSE, trim_ws = TRUE)
iran_equake = read_rds("data_eq/iran_earthquake.rds")

با استفاده از داده های زلزله ها در ایران و جهان به سوالات زیر پاسخ دهید.

۱. با استفاده از داده های historical_web_data_26112015.rds و استفاده از نمودار پراکنش سه بعدی بسته plotly نمودار طول، عرض و عمق زلزله ها را رسم نمایید. علاوه بر آن بزرگی هر نقطه را برابر بزرگی زمین لرزه قرار دهید. 

plot_ly(historical, x = ~Longitude, y = ~Latitude, z = ~Depth, size = ~Magnitude,
        marker = list( symbol = 'circle', sizemode = 'diameter'), sizes = c(1, 30))

۲. پویانمایی سونامی های تاریخی را بر حسب شدت بر روی نقشه زمین رسم نمایید.(از داده زلزله های بزرگ استفاده نمایید.) 

disaster %>%filter(FLAG_TSUNAMI == "Tsu" & !is.na(LATITUDE) & !is.na(LONGITUDE) & !is.na(LOCATION_NAME)& !is.na(EQ_PRIMARY)) %>% arrange(EQ_PRIMARY) -> t1
ggplot() + geom_polygon(data = map_data("world"), aes(long, lat, group=group), fill = "white", color = "lightblue") +
  geom_point(data = t1, aes(x = LONGITUDE,
                                     y = LATITUDE,
                                     frame = YEAR,
                                     color = EQ_PRIMARY,
                                     size = EQ_PRIMARY)) + 
  scale_color_continuous(low = "yellow", high = "red", guide = F) + 
  scale_size(guide = F) -> p
p

#gganimate(p, filename = "t.gif")

gganimate مشکل داشت و خروجی نمیداد!!!!

۳. نمودار چگالی دو بعدی زلزله های تاریخی ایران را رسم کنید.( از داده iran_earthquake.rds و لایه stat_density_2d استفاده نمایید). 

tehran_map = read_rds("data_eq/Tehrn_map_6.rds")
ggmap(tehran_map) + stat_density_2d(data = iran_equake , geom = "polygon", aes(x = Long , y = Lat , fill = ..level.. , alpha = ..level..)) + scale_alpha(range = c(0.4, 0.75), guide = FALSE) + scale_fill_gradient(low = "green", high = "red", guide = FALSE)

۴. احتمال اینکه در ایران در پنج سال آینده زلزله به بزرگی هفت ریشتر رخ دهد را محاسبه کنید. (از احتمال شرطی استفاده کنید.)

۵. بر اساس داده های زلزله های بزرگ ابتدا تعداد و متوسط کشته زلزله ها را بر حسب کشور استخراج نمایید. سپس نمودار گرمایی تعداد کشته ها را بر روی کره زمین رسم نمایید.(مانند مثال زیر!) 

disaster %>% filter(!is.na(TOTAL_DEATHS), !is.na(COUNTRY)) %>% group_by(COUNTRY) %>% summarise(TOTAL = sum(TOTAL_DEATHS), AVE = mean(TOTAL_DEATHS)) %>% mutate(COUNTRY_CODE = countrycode(COUNTRY , "country.name", "iso3c")) -> t4
hcmap(data = t4,joinBy = c("iso-a3", "COUNTRY_CODE"),value = "TOTAL",name = "TOTAL DEATHS")
hcmap(data = t4,joinBy = c("iso-a3", "COUNTRY_CODE"),value = "AVE",name = "AVERAGE DEATHS")

۶. با استفاده از داده لرزه های بزرگ و به وسیله طول، عرض، شدت، عمق مدلی برای پیش بینی تعداد کشته های زلزله بیابید. 

glm(TOTAL_DEATHS ~ LONGITUDE + LATITUDE + FOCAL_DEPTH + EQ_PRIMARY, disaster, family = "poisson") %>% summary()
## 
## Call:
## glm(formula = TOTAL_DEATHS ~ LONGITUDE + LATITUDE + FOCAL_DEPTH + 
##     EQ_PRIMARY, family = "poisson", data = disaster)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -384.38   -61.46   -39.20   -21.18  1521.82  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.389e+00  5.064e-03  -274.3   <2e-16 ***
## LONGITUDE    8.578e-04  6.032e-06   142.2   <2e-16 ***
## LATITUDE     1.973e-02  2.745e-05   718.9   <2e-16 ***
## FOCAL_DEPTH -2.631e-02  4.480e-05  -587.3   <2e-16 ***
## EQ_PRIMARY   1.348e+00  6.850e-04  1968.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 18100341  on 1051  degrees of freedom
## Residual deviance: 13160953  on 1047  degrees of freedom
##   (4974 observations deleted due to missingness)
## AIC: 13166049
## 
## Number of Fisher Scoring iterations: 8

۷. با استفاده از داده worldwide.csv به چند سوال زیر پاسخ دهید. تحقیق کنید آیا می توان از پیش لرزه، زلزله اصلی را پیش بینی کرد؟ 

worldwide$date = reshape2::colsplit(worldwide$time, pattern=" ", names=c("Part1", "Part2"))[,1]
worldwide %>% group_by(place, date) %>% mutate(mainmag = max(mag), ismain = (mag == max(mag))) %>% arrange(place,date,-mainmag) %>% mutate(timedif = as.numeric(first(time)-time)) %>% filter(timedif<0)->t7
sampleindex = sample(nrow(t7),nrow(t7)*0.9)
train = t7[sampleindex,]
test = t7[-sampleindex,]
fit1 = lm(data=train,mainmag~depth+mag)
predict1 = predict(fit1,test)
summary(fit1)
## 
## Call:
## lm(formula = mainmag ~ depth + mag, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.3577 -0.5962 -0.2049  0.4144  4.0499 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.176e+00  2.084e-02  104.41   <2e-16 ***
## depth       -8.814e-04  4.976e-05  -17.71   <2e-16 ***
## mag          6.730e-01  5.592e-03  120.35   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8615 on 32536 degrees of freedom
## Multiple R-squared:  0.3085, Adjusted R-squared:  0.3085 
## F-statistic:  7258 on 2 and 32536 DF,  p-value: < 2.2e-16
summary(predict1-test$mainmag)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -3.78070 -0.41130  0.19510 -0.00004  0.59146  1.35565
fit2 = lm(data=train,timedif~depth+mag)
predict2 = predict(fit2,test)
summary(fit2)
## 
## Call:
## lm(formula = timedif ~ depth + mag, data = train)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -55821 -18562   2446  20702  37117 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -43383.176    565.162 -76.762  < 2e-16 ***
## depth           -4.111      1.349  -3.047  0.00231 ** 
## mag           2580.180    151.633  17.016  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23360 on 32536 degrees of freedom
## Multiple R-squared:  0.008823,   Adjusted R-squared:  0.008762 
## F-statistic: 144.8 on 2 and 32536 DF,  p-value: < 2.2e-16
summary(predict1-test$timedif)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##     6.04 13964.08 32141.07 34521.58 53215.76 85840.84

۸. گزاره " آیا شدت زلزله به عمق آن بستگی دارد" را تحقیق کنید؟ (طبیعتا از آزمون فرض باید استفاده کنید.) 

worldwide %>% filter(!is.na(mag),!is.na(depth)) ->t8
cor.test(t8$depth,t8$mag)
## 
##  Pearson's product-moment correlation
## 
## data:  t8$depth and t8$mag
## t = 50.03, df = 60629, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1914584 0.2067469
## sample estimates:
##       cor 
## 0.1991147
ggplot(t8 %>% head(1000),aes(x = depth,y = mag)) + geom_point()

۹. میانگین سالانه زلزله ها را بر حسب کشور به دست آورید. آیا میتوان دلیلی در تایید یا رد تئوری هارپ ارائه کرد. 

worldwide %>% group_by(place,year) %>% summarise(avemag = mean(mag,na.rm = T))
## # A tibble: 925 x 3
## # Groups:   place [?]
##    place       year  avemag
##    <chr>       <chr>  <dbl>
##  1 Afghanistan 2015    4.33
##  2 Afghanistan 2016    4.35
##  3 Afghanistan 2017    4.30
##  4 Afghanistan 2018    4.35
##  5 Alabama     2016    2.75
##  6 Alabama     2017    2.66
##  7 Alaska      2015    2.99
##  8 Alaska      2016    3.08
##  9 Alaska      2017    3.11
## 10 Alaska      2018    3.18
## # ... with 915 more rows
ggplot() + geom_polygon(data = map_data("world"), aes(long, lat, group=group), fill = "white", color = "lightblue") +
  geom_point(data = worldwide %>% filter(mag > 4), aes(x = longitude,
                                     y = latitude,
                                     color = mag,
                                     alpha = .1)) + 
  scale_color_continuous(low = "yellow", high = "red", guide = F) + 
  scale_size(guide = F)

۱۰. سه حقیقت جالب در مورد زلزله بیابید.

آمار مرگ در اثر سونامی از زلزله بیشتر است. 

wilcox.test(
  disaster %>% filter(!is.na(FLAG_TSUNAMI) , !is.na(TOTAL_DEATHS)) %>% .$TOTAL_DEATHS,
  disaster %>% filter(is.na(FLAG_TSUNAMI) , !is.na(TOTAL_DEATHS)) %>% .$TOTAL_DEATHS,
  alternative = "greater"
)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  disaster %>% filter(!is.na(FLAG_TSUNAMI), !is.na(TOTAL_DEATHS)) %>%  and disaster %>% filter(is.na(FLAG_TSUNAMI), !is.na(TOTAL_DEATHS)) %>%     .$TOTAL_DEATHS and     .$TOTAL_DEATHS
## W = 317630, p-value = 8.633e-11
## alternative hypothesis: true location shift is greater than 0

مرگ و میر در ساعات کاری کمتر است. 

disaster %>% filter(!is.na(HOUR)) %>% group_by(HOUR) %>% summarise(count = n(),avemag = mean(EQ_PRIMARY,na.rm = T),avedeath = mean(TOTAL_DEATHS,na.rm = T)) -> t11
sum = sum(t11$count)
ggplot() + geom_bar(data = t11,aes(x = HOUR,y = avedeath),stat = "identity")